COMMENT ⊗   VALID 00009 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	.TITLE ARITH
C00004 00003	EXP & LOG -- FOR USE BY EXPONENTIATION ROUTINES & WORLD
C00011 00004	ASIN - FLOATING POINT, SINGLE PRECISION ARC-SINE FUNCTION
C00012 00005	ACOS - FLOATING POINT, SINGLE PRECISION ARC-COSINE FUNCTION
C00013 00006	ATAN2 - FLOATING POINT, SINGLE PRECISION ARC-TANGENT WITH TWO ARGUMENTS
C00016 00007	"SNCOS" - SINE/COSINE FUNCTION USING TABLE LOOKUP
C00019 00008	SINE/COSINE LOOK UP TABLES
C00023 00009	ARCTANGENT LOOK UP TABLES
C00027 ENDMK
C⊗;
.TITLE ARITH

;"SQRTF" - FLOATING POINT SQUARE ROOT

;COMPUTES THE SQUARE ROOT OF THE NUMBER LOADED INTO AC0 BY USING
;A LINEAR APPROXIMATION AND PERFORMING ONE CONVERGENCE INTERATION.
;THE ANSWER IS RETURNED IN AC0 AND IS ACCURATE TO APPROMIATELY
;+ OR - .0076%.   A SAMPLE CALLING SEQUENCE FOLLOWS:
;
;			LDF	NUMBER,AC0
;			JSR	PC,SQRTF
;
;EXECUTION TIME FOR THIS ROUTINE IS APPROXIMATELY 110 MICRO-SECONDS.

;REGISTERS USED:
;	AC0 PASSES ARGUMENTS AND AC1 IS GARBAGED
CODE
SQRTF:	ABSF	AC0		;ARG MUST BE > 0
	MOV	R0,-(SP)
	STEXP	AC0,R0		;EXPONENT ← EXPONENT/2
	STF	AC0,AC1
	INC	R0
	ASR	R0
	BCS	1$		
	LDEXP	#-1,AC1		;1/4 ≤ ARG < 1/2
	MULF	10$,AC1		;LINEAR APPROXIMATION
	ADDF	12$,AC1
	BR	2$
1$:	LDEXP	#0,AC1		;1/2 ≤ ARG ≤ 1
	MULF	11$,AC1
	ADDF	13$,AC1
2$:	LDEXP	R0,AC1		;SQRT ← ((ARG/GUESS)+GUESS)*.5
	DIVF	AC1,AC0
	MOV	(SP)+,R0
	ADDF	AC1,AC0
	MULF	#40000,AC0
	RTS	PC
	
DATA
10$:	.FLT2	0.8125000	;ASM
11$:	.FLT2	0.5781250	;ALG
12$:	.FLT2	0.3027340	;BSM
13$:	.FLT2	0.4218750	;BLG

;END OF "SQRTF"
;EXP & LOG -- FOR USE BY EXPONENTIATION ROUTINES & WORLD
;modified from those used by SAIL (GOGOL[S,AIL]/63P) by arg 9/78


;FLOATING POINT SINGLE PRECISION EXPONENTIAL FUNCTION
;THE ARGUMENT IS RESTRICTED TO THE FOLLOWING RANGE
;	-88.029693<X<88.029693
;IF X<-88.029693, THE PROGRAM RETURNS ZERO AS THE ANSWER
;IF X>88.029693, THE PROGRAM RETURNS ∞ AS THE ANSWER
;THE RANGE OF THE ARGUMENT IS REDUCED AS FOLLOWS:
;EXP(X) = 2**(X*LOG(E)BASE2) = 2**(M+F)
;WHERE M IS AN INTEGER AND F IS A FRACTION
;2**M IS CALCULATED BY ALGEBRAICALLY ADDING M TO THE EXPONENT
;OF THE RESULT OF 2**F. 2**F IS CALCULATED AS

;2**F = 2(0.5+F(A+B*F↑2 - F-C(F↑2 + D)**-1)**-1)

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
;	LDF X,AC0
;	JSR PC,EXP
;AC0 - AC3, R0 & R1 are used
;THE ANSWER IS RETURNED IN AC0

CODE
EXP:	STF AC0,AC1	;AC1 ← copy of argument
	ABSF AC1	;AC1 ← absolute value of argument
	CMPF 10$,AC1	;Make sure argument is in proper range
	CFCC
	BGE 2$		;Skip if it's okay
	TSTF AC0	;See which end of the range it exceeds
	CFCC
	BMI 1$
	LDF INF,AC0	;Too big - use largest number possible
	RTS PC
1$:	CLRF AC0	;Too small - return zero
	RTS PC

2$:	MODF 15$,AC0	;Multiply arg by log e (base 2) - AC0 ← fractional part
			;				  AC1 ← integer part
	STCFI AC1,R0	;R0 ← M
	STF AC0,AC1	;Copy F
	MULF AC1,AC1	;Form F↑2
	LDF 12$,AC2	;Get B
	MULF AC1,AC2	;AC2 ← B * F↑2
	ADDF 14$,AC1	;AC1 ← F↑2 + D
	LDF 13$,AC3	;AC3 ← C
	DIVF AC1,AC3	;AC3 ← C/(F↑2 + D)
	SUBF AC3,AC2	;AC2 ← B * F↑2 - C/(F↑2 + D)
	SUBF AC0,AC2	;AC2 ← B * F↑2 - F - C/(F↑2 + D)
	ADDF 11$,AC2	;AC2 ← A + B * F↑2 - F - C/(F↑2 + D)
	DIVF AC2,AC0	;AC0 ← F / (A + B * F↑2 - F - C/(F↑2 + D)) = 2 ↑ (F-1)
	ADDF HALF,AC0	;ADD 0.5
	STEXP AC0,R1	;R1 ← exponent of 2 ↑ (F-1)
	INC R1		;R1 ← exponent of 2 ↑ F
	ADD R0,R1	;Add in M - i.e. 2↑M * 2↑F
	LDEXP R1,AC0	;Load correct exponent
	RTS PC		;& return
DATA
10$:	.FLT2	88.029693	;EXP(88.029693) = INF
11$:	.FLT2	9.95459578
12$:	.FLT2	0.03465735903
13$:	.FLT2	617.97226953
14$:	.FLT2	87.417497202
15$:	.FLT2	1.44269504088	;LOG(E), BASE 2
HALF::	.FLT2	0.5
INF::	.WORD	077777, 177777	;Largest possible positive floating point number

;FLOATING POINT SINGLE PRECISION LOGARITHM FUNCTION
;LOG(ABSF(X)) IS CALCULATED BY THE SUBROUTINE, AND AN
;ARGUMENT OF ZERO IS RETURNED AS MINUS INFINITY. THE ALGORITHM IS

;LOGE(X) = (I + LOG2(F))*LOGE(2)
;WHERE X = F * 2↑I, AND LOG2(F) IS GIVEN BY
;LOG2(F) = C1*Z + C3*Z↑3 + C5*Z↑5 - 1/2
;AND Z = (F-SQRT(2))/(F+SQRT(2))

;THE CALLING SEQUENCE FOR THE ROUTINE IS AS FOLLOWS:
;	LDF X,AC0
;	JSR PC,EXP
;AC0 - AC2 & R0 are used
;THE ANSWER IS RETURNED IN AC0
CODE
LOG:	ABSF AC0	;Make sure argument is >0
	CFCC
	BNE 1$
	LDF INF,AC0	;If arg is 0 return largest possible negative number
	NEGF AC0
	RTS PC
1$:	CMPF 10$,AC0	;Check for 1.0 argument
	CFCC
	BNE 2$
	CLRF AC0	;If arg is 1.0 return 0
	RTS PC

2$:	STEXP AC0,R0	;R0 ← Exponent of argument
	LDEXP #1,AC0	;AC0 ← F
	STF AC0,AC1	;Copy F
	SUBF 11$,AC0	;AC0 ← F - Sqrt(2)
	ADDF 11$,AC1	;AC1 ← F + Sqrt(2)
	DIVF AC1,AC0	;AC0 ← (F-Sqrt(2))/(F+Sqrt(2)) = Z
	STF AC0,AC1	;Copy Z
	MULF AC1,AC1	;AC1 ← Z↑2
	LDF 12$,AC2	;AC2 ← C5
	MULF AC1,AC2	;AC2 ← C5 * Z↑2
	ADDF 13$,AC2	;AC2 ← C3 + C5 * Z↑2
	MULF AC1,AC2	;AC2 ← C3 * Z↑2 + C5 * Z↑4
	ADDF 14$,AC2	;AC2 ← C1 + C3 * Z↑2 + C5 * Z↑4
	MULF AC0,AC2	;AC2 ← C1 * Z + C3 * Z↑3 + C5 * Z↑5 = LOG2(F) + 1/2
	SUBF HALF,AC2	;AC2 ← LOG2(F)
	LDCIF R0,AC0	;AC0 ← I
	ADDF AC2,AC0	;AC0 ← LOG2(X)
	MULF 15$,AC0	;AC0 ← LOG2(X) * LN(2) = LN(X)
	RTS PC		;Done - return

;CONSTANTS
DATA
10$:	.FLT2 1.0
11$:	.FLT2 1.414213562374	;Square root of 2
12$:	.FLT2 0.5989786496	;C5
13$:	.FLT2 0.9614706323	;C3
14$:	.FLT2 2.8853912903	;C1
15$:	.FLT2 0.69314718056	;Natural log of 2

;ASIN - FLOATING POINT, SINGLE PRECISION ARC-SINE FUNCTION

;COMPUTES THE ARCSINE OF AN VALUE USING THE FOLLOWING ALGORITHM
;
;	ASIN(X) = ATAN2(X,SQRT(1-X**2))
;
;WHERE X IS IN THE RANGE -1.0 TO +1.0.  THE ARGUMENT X MUST BE 
;LOADED INTO AC0 BEFORE THE CALL TO ASIN.  AFTER EXECUTION ASIN
;RETURNS THE ARCSINE VALUE IN AC0 ( IN DEGREES ).

;REGISTERS USED:
;
;	AC0 PASSES ARGUMENTS
;	AC1,AC2  GARBAGED
CODE
ASIN: 	STF	AC0,AC2		;COMPUTE SQRT(1-SIN↑2)
	MULF	AC0,AC0
	NEGF	AC0
	ADDF	#40200,AC0
	JSR	PC,SQRTF
	STF	AC0,AC1
	STF	AC2,AC0
	JSR	PC,ATAN2 	;ARC TANGENT( SIN, SQRT(1-SIN↑2) )
	RTS	PC

;END OF "ASIN"
;ACOS - FLOATING POINT, SINGLE PRECISION ARC-COSINE FUNCTION

;COMPUTES THE ARC-COSINE OF AN VALUE USING THE FOLLOWING ALGORITHM
;
;	ACOS(X) = ATAN2( SQRT(1-COSX↑2),COS X )
;
;WHERE X IS IN THE RANGE -1.0 TO +1.0.  THE ARGUMENT X MUST BE 
;LOADED INTO AC0 BEFORE THE CALL TO ACOS.  AFTER EXECUTION ACOS
;RETURNS THE ARC-COSINE VALUE IN AC0 ( IN DEGREES ).

;REGISTERS USED:
;
;	AC0 PASSES ARGUMENTS
;	AC1,AC2  GARBAGED

CODE
ACOS: 	STF	AC0,AC2		;COMPUTE SQRT(1-COS↑2)
	MULF	AC0,AC0
	NEGF	AC0
	ADDF	#40200,AC0
	JSR	PC,SQRTF
	STF	AC2,AC1
	JSR	PC,ATAN2 	;ARC TANGENT( SQRT(1-COS↑2), COS )
	RTS	PC

;END OF "ACOS"
;ATAN2 - FLOATING POINT, SINGLE PRECISION ARC-TANGENT WITH TWO ARGUMENTS

;COMPUTES THE ARC-TANGENT OF A/B USING A TABLE LOOK UP SCHEME.  SINCE
;TWO ARGUMENTS ARE USED SINGULARITIES ARE AVOIDED AT MULTIPLES OF PI/2 
;AND THERE IS NO AMBIGUITY CONCERNING QUADRANTS.  THE ARGUMENT A MUST
;BE LOADED INTO AC0 AND B INTO AC1 BEFORE CALLING ATAN2.  AFTER EXECUTION,
;ATAN RETURNS THE ARC-TANGENT IN DEGREES IN AC0.  RUN TIME IS APPROXIMATELY
;140 MICRO SECONDS.

;REGISTERS USED:
;
;	AC0,AC1 PASSES ARGUMENTS
;	AC2 GARBAGED

SSIN==1		;SINE IS NEGATIVE
SCOS==2		;COSINE IS NEGATIVE
CMPANG==4	;COMPLEMENTARY ANGLE

ATAN2: 	MOV	R0,-(SP)	;SAVE REGISTERS
	MOV	R1,-(SP)
	CLR	R0
	TSTF	AC0		;SIN > 0 ?
	CFCC
	BLT	SINNEG
	BGT	SINPOS
	TSTF	AC1		;SIN = 0, COS > 0?
	CFCC
	BPL	.+6
    	LDF	#42064,AC0	;THETA = 180 DEGREES
JRET:	MOV	(SP)+,R1
	MOV	(SP)+,R0
	RTS	PC

SINNEG:	NEGF	AC0		;ABS(SIN)
	BIS	#SSIN,R0  	;INDICATE SIN < 0
SINPOS:	TSTF	AC1		;COS > 0 ?
	CFCC
	BLT	COSNEG
	BGT	COSPOS
	LDF	#41664,AC0	;THETA = 90 DEGREES
	BIT	#SSIN,R0
	BEQ	JRET
	NEGF	AC0
	BR	JRET

COSNEG:	NEGF	AC1		;ABS(COS)
	BIS	#SCOS,R0 	;INDICATE COS < 0
COSPOS:	CMPF	AC0,AC1
	CFCC
	BGT	1$		;MUST EXCHANGE
	BNE	2$
	LDF	#41464,AC0	;THETA = 45 DEGREES
	BR	4$

1$:	DIVF	AC0,AC1
	BIS	#CMPANG,R0
	STF	AC1,AC0
	BR	3$

2$:	DIVF	AC1,AC0
3$:	MODF	#42000,AC0	;GET INDEX, INTERPOLATION FACTOR
	STCFI	AC1,R1		;INDEX
	ASH	#2,R1
	LDF	ARCTAN+4(R1),AC2	;INTERPOLATE
	LDF	ARCTAN(R1),AC1
	SUBF	AC1,AC2
	MULF	AC0,AC2
	STF	AC1,AC0
	ADDF	AC2,AC0

4$:	BIT	#CMPANG,R0	;COMPLEMENT?
	BEQ	5$
	SUBF	#41664,AC0
	NEGF	AC0
5$:	BIT	#SCOS,R0	;COS THETA < 0?
	BEQ	6$
	SUBF	#42064,AC0
	NEGF	AC0
6$:	BIT	#SSIN,R0	;SIN THETA < 0?
	BEQ	JRET
	NEGF	R0
	BR	JRET

;END OF "ATAN2"
;"SNCOS" - SINE/COSINE FUNCTION USING TABLE LOOKUP

;THIS PROGRAM CALCULATES BOTH THE SINE AND THE COSINE OF A ANGLE USING
;A TABLE LOOP UP PROCEDURE.  THE IMPLEMENTED APPROXIMATION EQUATIONS
;ARE AS FOLLOWS:
;		SIN(X) = SIN(A) + (B/I)*[SIN(A+I)-SIN(A)]
;		COS(X) = COS(A) + (B/I)*[COS(A+I)-COS(A)]
;	WHERE
;		I = 90/128 DEGREES
;		A = INTEGER(X*128/90)
;		B = REMAINDER(X*128/90)
;
;THE ANGLE SHOULD MUST BE IN DEGREES AND MUST BE LOADED INTO AC0 BEFORE
;CALLING THIS ROUTINE.  ON EXITING, THE SIN IS LEFT IN AC0 AND THE COSINE
;IS LEFT IN AC1.  EXECUTION TIME IS APPROXIMATELY 120 MICRO-SECONDS.

;REGISTERS USED:
;	AC0,AC1 PASS ARGUMENTS AND ARE ALTERED
;	AC2,AC3 ARE GARBAGED

SNCOS:	MODF	3$,AC0		;SEPARATE INDEX AND FRACTION
	MOV	R1,-(SP)	;SAVE REGISTERS
      	MOV	R0,-(SP)	
	STCFI	AC1,R0		;INDEX INTO LOOKUP TABLE
	MOV	R0,R1		;SAVE SIGN OF ANGLE
	BIC	#177600,R0
	ASH	#2,R0
	LDF	SINTAB+4(R0),AC1	;SIN(A+I)
	LDF	SINTAB(R0),AC2	;SIN(A)
	SUBF	AC2,AC1
	ASR	R1
	MULF	AC0,AC1		;(B/I)*[SIN(A+1)-SIN(A)]
	SWAB	R1
	ADDF	AC2,AC1		;NOW HAVE ABS(SIN(X))
	NEG	R0		;ENTER TABLE FROM OPPOSITE END
	LDF	COSTAB-4(R0),AC2	;COS(A+I)
	LDF	COSTAB(R0),AC3	;COS(A)
	SUBF	AC3,AC2
	MULF	AC0,AC2		;(B/I)*[COS(A+I)-COS(A)]
	MOV	(SP)+,R0	;DON'T NEED THIS ANY MORE
	BIT	#40000,R1	;WHICH IS THE SIN?
	ADDF	AC3,AC2		;NOW HAVE ABS(COS(X))
	BEQ	1$
	STF	AC2,AC0		;SWITCH SIN ↔ COSINE
	BR	2$
1$:	STF	AC1,AC0		;SAVE SIN, COSINE
	STF	AC2,AC1
2$:	TST	R1    		;SIN ← -SIN IF X IN QUAD 3 OR 4
	BPL	.+4
	NEGF	AC0
	ADD	#40000,R1   	;COS ← -COS IF QUADRANT 2 OR 3
	BPL	.+4
	NEGF	AC1
	MOV	(SP)+,R1
	RTS	PC
DATA
3$:	.FLT2	1.422222	;128/90 = 1.422222    

;END OF "SNCOS"
;SINE/COSINE LOOK UP TABLES
DATA
SINTAB:	.WORD	     0,     0, 36511,  7220, 36711,  5260, 37026,141454
	.WORD	 37110,175460, 37173, 25564, 37226,124405, 37257,133200
	.WORD	 37310,136466, 37341,136056, 37372,131163, 37411,147606
	.WORD	 37426, 40203, 37442,125666, 37457, 10242, 37473, 67317
	.WORD	 37507,142702, 37524, 12401, 37540, 56023, 37554,115177
	.WORD	 37570,147714, 37602, 76700, 37610,107223, 37616,115042
	.WORD	 37624,120062, 37632,120206, 37640,115345, 37646,107422
	.WORD	 37654, 76324, 37662, 61757, 37670, 42052, 37676, 16512
	.WORD	 37703,167425, 37711,134523, 37717, 75712, 37725, 33101
	.WORD	 37732,164201, 37740,111117, 37746, 31565, 37753,145673
	.WORD	 37761, 55352, 37766,160313, 37774, 56447, 40000,163744
	.WORD	 40003,116075, 40006, 45603, 40010,172633, 40013,115153
	.WORD	 40016, 34732, 40020,151715, 40023, 64052, 40025,173331
	.WORD	 40030, 77700, 40033,  1306, 40035, 77721, 40037,173313
	.WORD	 40042, 63631, 40044,151045, 40047, 33126, 40051,112025
	.WORD	 40053,165512, 40056, 35736, 40060,102673, 40062,144311
	.WORD	 40065,  2363, 40067, 35043, 40071, 64102, 40073,107473
	.WORD	 40075,127371, 40077,143547, 40101,154160, 40103,161001
	.WORD	 40105,162003, 40107,157145, 40111,150422, 40113,135770
	.WORD	 40115,117402, 40117, 75037, 40121, 46475, 40123, 14111
	.WORD	 40124,155461, 40126,112745, 40130, 44123, 40131,171152
	.WORD	 40133,112032, 40135, 26523, 40136,137005, 40140, 43041
	.WORD	 40141,142630, 40143, 36132, 40144,125131, 40146,  7610
	.WORD	 40147, 65730, 40150,137474, 40152,  4647, 40153, 45414
	.WORD	 40154,101537, 40155,131223, 40156,154236, 40157,172563
	.WORD	 40161,  4410, 40162, 11522, 40163, 12110, 40164,  5735
	.WORD	 40164,175013, 40165,157306, 40166,135007, 40167,105705
	.WORD	 40170, 51770, 40171, 11230, 40171,143635, 40172, 71402
	.WORD	 40173, 12276, 40173,126315, 40174, 35450, 40174,137711
	.WORD	 40175, 35254, 40175,125713, 40176, 11443, 40176, 70260
	.WORD	 40176,142155, 40177,  7130, 40177, 47155, 40177,102253
	.WORD	 40177,130417, 40177,151627, 40177,166103, 40177,175421
COSTAB:	.WORD	 40200,     0
;ARCTANGENT LOOK UP TABLES

ARCTAN:	.WORD	     0,     0, 37745, 26657, 40145, 25033, 40253,155433
	.WORD	 40345, 15712, 40417, 25252, 40453,141371, 40510, 52757
	.WORD	 40544,161252, 40600,132045, 40616,171370, 40635, 26536
	.WORD	 40653, 61353, 40671,111461, 40707,136703, 40725,161067
	.WORD	 40744,    42, 40762, 13434, 41000, 11537, 41007, 13517
	.WORD	 41016, 13456, 41025, 11311, 41034,  4757, 41042,176160
	.WORD	 41051,165034, 41060,151305, 41067,133074, 41076,112127
	.WORD	 41105, 66347, 41114, 37702, 41123,  6276, 41131,151661
	.WORD	 41140,112164, 41147, 47337, 41156,  1313, 41164,130024
	.WORD	 41173, 53224, 41200,175425, 41204, 43531, 41207,110006
	.WORD	 41212,152413, 41216, 13131, 41221, 51741, 41224,106624
	.WORD	 41227,141545, 41232,172506, 41236, 21451, 41241, 46403
	.WORD	 41244, 71310, 41247,112153, 41252,130743, 41255,145445
	.WORD	 41260,160046, 41263,170336, 41266,176504, 41272,  2517
	.WORD	 41275,  4367, 41300,  4065, 41303,  1402, 41305,174527
	.WORD	 41310,165457, 41313,154203, 41316,140516, 41321,122614
	.WORD	 41324,102472, 41327, 60121, 41332, 33320, 41335,  4263
	.WORD	 41337,152767, 41342,117232, 41345, 61232, 41350, 20767
	.WORD	 41352,156256, 41355,111276, 41360, 42050, 41362,170354
	.WORD	 41365,114410, 41370, 36176, 41372,155520, 41375, 72574
	.WORD	 41400,  2603, 41401, 46770, 41402,112035, 41403,153763
	.WORD	 41405, 14574, 41406, 54272, 41407,112655, 41410,150127
	.WORD	 41412,  4272, 41413, 37330, 41414, 71263, 41415,122115
	.WORD	 41416,151647, 41420,   304, 41421, 25647, 41422, 52122
	.WORD	 41423, 75306, 41424,117410, 41425,140431, 41426,160374
	.WORD	 41427,177264, 41431, 15103, 41432, 31655, 41433, 45365
	.WORD	 41434, 60034, 41435, 71450, 41436,102032, 41437,111366
	.WORD	 41440,117676, 41441,125167, 41442,131443, 41443,134706
	.WORD	 41444,137143, 41445,140375, 41446,140631, 41447,140070
	.WORD	 41450,136337, 41451,133621, 41452,130121, 41453,123444
	.WORD	 41454,116014, 41455,107415, 41456,100053, 41457, 67551
	.WORD	 41460, 56313, 41461, 44125, 41462, 31011, 41463, 14753
	.WORD	 41464,     0


;END OF "ARITH"